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Abstract 

We integrate the one-dimensional nonlinear Schrodinger equation numerically for 
solitons moving in external potentials. In particular, we study the scattering off an 
interface separating two regions of constant potential modeled by a linear ramp. 
Transmission coefficients and inelasticities are computed as functions of the poten- 
tial difference and the slope of the ramp. Our data show that the ramp's slope 
has a strong influence revealing unexpected windows of reflection in a transmis- 
sion regime. The transmission coefficients for very small potential differences are 
compared with the theoretical predictions derived by perturbation theory. Also the 
time evolution of the solitary waves after the scattering is studied. We observed 
that they in general behave like solitons with an amplified amplitude. Due to this, 
they oscillate. The oscillation period is measured and compared with theoretical 
predictions. 

PACS number(s): 03.40.Kf, 0.3.65. Ge, 42.65.Tg, 42.81. Dp 
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1 Introduction 

Recent years have seen a considerable growth in the interest for nonlinear 
partial differential equations with soliton solutions. The most important and 
very well studied equations are the Korteweg-de Vries equation with serveral 
of its modifications, the sine-Gordon equation and the nonlinear Schrodinger 
equation (NLSE). Historically, the first person who observed a solitary wave 
was the Scotsman J.S. Russell [1,2]. He discovered on the Edinburgh- Glasgow 
channel that a large solitary 'heap of water', as he calls this phenomenon, 
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moves over a very long distance without been dispersed. The equation of mo- 
tion to describe these waves in shallow water was first found by D.J. Korteweg 
and G. de Vries in 1895 [3]. After this it took some time, until 1967, when 
C.S. Gardener, J.M. Greene, M.D. Kruskal and R.M. Miura [4,5] showed that 
the Korteweg-de Vries equation can be solved analytically with the so-called 
inverse scattering theory (1ST). Moreover, it was also possible to solve the 
sine-Gordon equation and the NLSE with the aid of the same method [6-8]. 

In the present paper we study the NLSE which appears with its variants in 
problems drawn from disciplines as diverse as optics, solid state, particle and 
plasma physics. There, the NLSE describes phenomena such as modulational 
instability of deep water waves [9] , propagation of heat pulses in anharmonic 
crystals, helical motion of very thin vortex filaments, models of protein dynam- 
ics [10], nonlinear modulation of collisionless plasma waves [11], transmission 
of pulses through junctions in optical fibers [12-14], and self-trapping of light 
beams in optically nonlinear media [15,16]. In the last decade the optical com- 
munication via solitons has become a field of special interest [17-20]. In all 
these applications, the main interest is due to the fact that the NLSE is the 
most fundamental model which describes the interplay between weak disper- 
sion and nonlinearity for a wave envelope. The existence of soliton solutions is 
one consequence of this. These are solitary waves with well defined pulse-like 
shapes and remarkable stability properties [21]. 

A great deal of current interest is directed to the question how solitons be- 
have under the influence of external perturbations [15,22-27]. The perturba- 
tions describe e.g. various dissipative effects, inhomogeneities of a medium 
or effects of external fields. They can be naturally divided into two classes: 
Hamiltonian and dissipative. Unfortunately, most of them destroy the inte- 
grability of the NLSE. But in many cases, the perturbations are small and 
the techniques of analytical perturbations theory can be applied to extract 
some results. Perturbation-induced effects are of interest mainly because they 
represent physical phenomena that cannot be comprised by exactly integrable 
models. 

Here we shall limit ourselves to such perturbations which can be described 
by potentials. They preserve the Hamiltonian structure of the NLSE [11], 
but not its complete integrability. More precisely, we shall consider the one- 
dimensional NLSE with external potentials which are constant outside a finite 
interval. Thereby we emulate the effect of an interface with finite width be- 
tween two different media in which the solitons have different characteristics. 
We study initial conditions consisting of one single soliton which moves to- 
wards the interface. There it will interact and, in general, it will not just be 
transmitted or reflected. Instead, we expect that the soliton will also undergo 
inelastic scatterings where it breaks up either into several solitons or into 
nonsolitary waves, or both. 
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This problem has been studied previously by several authors. While pertur- 
bative approaches were used in [15,16,25,28], straightforward numerical inte- 
grations were made in [26,29,27]. Both approaches showed that the soliton 
behaves just like a Newtonian particle if the force created by the potential is 
sufficiently weak. This is to be expected, but the problem what happens when 
the force is strong was not studied extensively. The problem with discontin- 
uous potentials is considered in [30,31], where the main interest was broken 
conservation laws and soliton tunneling. In [32] the authors have studied the 
behavior of the soliton scattering at a potential ramp only qualitatively. The 
main goal was the successful application of higher order symplectic integrators 
to the NLSE. The spatial inhomogeneity considered in [16] was more compli- 
cated, because the authors allowed also a discontinuity in the nonlinear part. 
That lead to a different behavior of the scattered solitons. 

It is one purpose of the present paper to study the potential ramp [28,33] 
in a more quantitative and systematical way. The main reason to investigate 
this kind of external potential was the fact that it reveals a very interesting 
feature. Namely, while for most parameter values of the slope and the hight 
of the potential the soliton behaves as expected in a "particle-like" manner, 
for special sets of parameters the soliton shows decidedly inelastic proper- 
ties. In particular, the nonlinear reflection and transmission coefficients of the 
mass and of the energy are shown to possess windows of behavior that are 
quite counter-intuitive at first sight. We find that as the width of the poten- 
tial ramp is varied, the transmission coefficient abruptly drops from a finite 
non-zero value to near zero and then just as abruptly back to a value near 
1. This feature does not disappear if one smoothes the potential. Only for 
very smooth potentials the behavior of the soliton equals more and more that 
of a Newtonian particle in an equivalent potential. First we have compared 
some of our numerical results for the flat potential step with predictions by 
perturbation theory. They agree very well in the expected validation range. 

We also have extensively discussed the influence on the soliton parameters 
(amplitude and velocity) during the scattering process. There we come to 
the conclusion that the solitons escaping in both regions after the scattering 
are "amplified solitons" . These are solitons whose amplitude is multiplied by 
some constant factor. Their behavior was studied by Kath and Smyth [34] 
who derived a formula to predict the amplitude oscillation period depending 
on the soliton parameters. This formula is applied to our data and is in quite 
good agreement with them. 

This article is organized as follows. In the next section we briefly summarize 
some basic relations and equations concerning the soliton solutions of the 
NLSE and introduce our problem. A detailed discussion of the transmission 
coefficients is made in Sec. 3. There some simulations are compared with 
results from the perturbation analysis, as far as this is possible within its 
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validity range. In Sec. 4 we take a closer look at the shape of the soliton after 
it is scattered at the inhomogeneity. An amplified soliton ansatz is made and 
the parameters of the corresponding solitary wave are determined. A summary 
and outlook is given in Sec. 5. In the appendix the applied integration method 
(symplectic integrators) is shortly described. The advantage over the other 
explicit integration schemes commonly used in this field, such as Runge-Kutta 
or predictor-corrector methods, is explained. 



2 The NLSE soliton solution 



Using appropriate units, we can write the NLSE as 

.dW(x,t) ld 2 ty(x,t) |T . , i9t . , Tr , N T . . 
,_iJ^ = --—Aj^ - \^(x,t)\ 2 ^(x,t) + V(x)^(x,t), (1) 



where V(x) is the external potential. We shall use for the latter a piecewise 
linear ansatz, with V(x) = for x < 0, V(x) = V > for x > x > 0, and 
linearly rising for x between and xq, 



V(x) 





xVq/xq 

V 



x < 

< X < Xq 
X > Xq 



(2) 



We call the negative x-axis region /, while region II is the region x > (where 
V(x) > 0). 

The NLSE with the perturbation (2) is not completely intergrable. But for 
a constant potential Vo it is, and the soliton solutions of Eq.(l) form a two- 
parameter manifold (apart from translations). Taking as parameters the ve- 
locity v and the amplitude a, these solutions read [35] 

( x> t ) = - . " e i{v*+U° a -*)/*-v ]t} (3) 



We denote the velocity of the incoming soliton as v - Using a suitable rescaling 
of x, t and ^, we can always choose its amplitude as ao = 1/2, without loss of 
generality. 

Among the infinitely many conserved quantities (for V(x) = const!) [25], the 
following three are of particular interest: 
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(i) the mass or "number of quasi-particles" 

N — J |^| 2 dx , 

which is conserved due to phase invariance; 

(ii) the energy 



(4) 



E 



dx 



- \\^\ 4 + V(x) l^l 2 ) <!•'■ 



(5) 



which is conserved due to time invariance; 
(iii) and the momentum 



P 1 f (V** 9 * \p d ^*\ d 
2i J \ dx dx 



(6) 



which is conserved due to translation invariance. 



For the soliton given by Eq.(3) this yields 

N = 2a, (7a) 
P = vN and (7b) 
E = (v 2 /2-a 2 /Q)N + (V)N , (7c) 

where the average over V(x) is taken with weight oc \^/\ 2 as indicated by 
Eq.(5). 

In the cases where the soliton retains its shape, the motion of the soliton can 
be understood as that of a Newtonian particle moving in the presence of an 
effective potential [36]. This assumption is valid when the force created by 
the external potential is sufficiently weak. In our case this is true when the 
local disturbance is small, i.e. Vq <C K or when x ^ i n our case. Then 
we can apply Ehrenfest's theorem to get the effective potential describing the 
motion of the Newtonian particle with mass m = 2a = 1 and initial velocity 
vq. Thus, in the case where the potential is much smaller than the kinetic 
energy of the particle (Vo <C K ,K = v 2 /2) we expect transmission with a 

final velocity i>zr i0 ut = \J v o ~ 2Vo- This follows from the energy conservation 
law and the assumption that the soliton remains entire and no radiation is 
emitted. This also predicts that there is no transmission if vq < y/2V (i.e., 
Kq < Vq) and x ^ V . But we will see that this is more or less the exception 
in our simulations and the true behavior will lead to unexpected windows in 
the transmission coefficients. 

To get some more insight into the problem, let us first look at the motion of 
a Newtonian particle with a shape \*tf(x, 0)| 2 in a potential (2). We define the 
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centre of mass of the extended particle 

1 °° 

z{t) = N f dxx ^M\ 2 ( 8 ) 

—00 



and derive an effective potential according to which the centre of mass is 
moving. Starting from the equation of motion we get: 



—00 
1 x ° v 

= -- fdx^(x,t)\\ (9) 
2a J x Q 




where we have chosen the potential from Eq. (2). While this is still exact 
for solitary waves, we now approximate the wave function by an unperturbed 
soliton centered at z(t). The shape of the soliton is then given by 



a 2 



l*(M)r = — T3H TTyTT ■ (10) 

cosh [a(x — z(t))\ 



We are aware that this assumption is valid only if the local derivative of the 
potential is very small. 



From the Eqs.(9) and (10) we can derive the effective potential 

cosh(az) 



V eS (z) 



Vo 

Xq 



In 



s cosh(a(,2 — Xq)) t 



Xq 

2 



;n) 



according to which the centre of mass of a soliton (3) is moving. In the limit 
clxq ^> 1 one sees that V e g(z) — > V(z). In the upper limit xq — * in which the 
potential becomes a step function the effective potential takes the form 

V eS {z) = ^(tanh(az) + 1) , (12) 



which agrees with the potential derived by Nogami and Toyama [29]. 

As one can easily see these functions are smooth sigmoids in contrast to V(x). 
Due to this it is possible that the soliton's centre of mass penetrates into the 
region II (i.e., where z > xq) although Vq > K . 
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Because of the fact that the effective potential has no local minima the soliton 
cannot be at rest at any place (provided that ?; ^ 0). This statement does not 
hold if one also allows a change in the nonlinear term at the inhomogeneity 
[15]. Then it is possible that the soliton is trapped in the vicinity of the 
inhomogeneity. 



3 Transmission and reflection coefficients 

Within our simulations we have intensively studied the transmission coef- 
ficients of the mass, the momentum, and the energy. We found that they 
strongly depend on the geometry of the external potential. Especially the 
slope of the ramp has a very strong influence on them, which is not expected 
for the equivalent particle method [36]. For the integration we used a fourth 
order symplectic integration algorithm [32] which is described in the appendix. 

We measured the mass Ni and the energy Ei in each region (i = I, II) sepa- 
rately. Notice that iV and E are conserved even in the presence of a noncon- 
stant potential [30]. Thus, they are of special importance among the infinite 
numbers of conserved quantities in the case where the potential is constant 
[25] . Accordingly, we can define two sets of transmission and reflection coeffi- 
cients. We call them T/v, Rn and Tp, Re, 

T N = , R N = -jjr = 1 - T N (13) 

and 

T E = ^T' R E = jr = l-T E . (14) 

A physical interpretation of T N can be the following: if A" is the number 
of quasi-particles of the incoming soliton, then is T/v the quantity of quasi- 
particles penetrating the potential; resp. is T E the energy of these quasi- 
particles. 

Another quantity of interest is the momentum. It is only conserved for a 
constant potential; otherwise the translation invariance is broken. Nevertheless 
the momentum of the transmitted and reflected wave can be defined in the 
same way 

Tp = ^ , Rp = § • (15) 
Due to nonconservation, in general we have Tp + Rp ^ 1. 
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All these coefficients are functions depending on the incident velocity vq of 
the soliton, the height Vq and the length x$ of the potential ramp. Our main 
interest lays in the detailed dependence on these parameters. Thus, we numer- 
ically evaluate the coefficients for a few initial velocities of the soliton varying 
the length and the height over a wide range. 

Primarily, we will study these coefficients in a parameter range where the 
perturbation approach is not valid. Only very few works in this field have 
dealt with perturbations which are of the same order as the energy of the 
soliton. Up to now the only way to get some insight is to simulate the NLSE 
numerically. 

In the following simulations we used At = 0.005, Ax = ^ ~ 0.2, and a = 1/2. 
Hence the initial soliton takes the form 

V(x,t = 0) = {2cosh[i(a;-x off )]}~ 1 e^ o:r , (16) 

where the offset x Q g is the position of the soliton at t — 0. The kinetic energy 
of the soliton is, according to Eq. (5), K = Vq/2. 

In addition to the transmission and reflection coefficients we also registered all 
local maxima of |\l/(a;)| 2 with |\I/(:r)| 2 > 1/3000. The threshold was chosen in 
such a way that no essential information is lost, but small radiation and the 
noise which is produced due to the numerical integration are not taken into 
account. 



3.1 Transmitted mass and energy 



3.1.1 Potential step 

Let us first examine the potential step (xo = 0), because there it is possible to 
compare the numerically measured transmission coefficients with the results 
calculated by perturbation theory. These calculations are valid only in two 
situations: when the potential step is very low in comparison to the kinetic 
energy of the soliton, and when the potential is much higher than the kinetic 
energy (reflection at a hard wall). 

In the first case the soliton is transmitted nearly without loss of radiation. The 
determination of the transmission coefficient is based on the 1ST technique 
using a first order Born approximation. I.e., the velocity of the soliton is 
assumed to be constant during the scattering process. This is equivalent to 
the inequality [24] 

K > 8a 2 \V \ . (17) 
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In order to calculate analytically the transmission coefficient for a single soliton 
of the perturbed NLSE (V < K ) 



.dV(x,t) ld 2 V{x,t) 



dt 



2 dx 2 



\^(x,t)\ 2 ^(x,t) + V f(x)^(x,t), (18) 



we follow the derivation given in [37]. There the authors study a soliton which 
is scattered at an impurity (/ is a Dirac (^-distribution), whereas in our case 
the function / is the step function 



9(x) 



: x < 

1 : x > 



(19) 



A calculation with this potential analogous to that in [37] gives a formula for 
the amount of radiation which is emitted in the backward direction, 



u 

R N = Jd\n rad (X) , 



(20) 



and the energy of that radiation, 



R E = /dA4A 2 n rad (A), 



(21) 



Here n ra d(A) denotes the total density of linear waves emitted by a soliton 
during the scattering, 



irV 2 

n rad (A) = ^f (3 2 (\) sech 2 



7T 



a v 



2X Z + — - 

CLV \ 2 



oo oo 



x Jdx J dx' f(x) f(x') exp[iP(\)(x - x')} , 



-oo — oo 



where 



m = - 



v \ 2 a 2 



A) 4 



The analytic calculation yields for the step function 



R N = ^4- Id A sech 2 



7r / x9 a 2 w 2 " 
— 2A 2 + — - — 
aw 1 2 8 



(22) 



(23) 
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Fig. 1. The transmission coefficients and Tg vs. Vq/Kq with fixed xq = and 
vq = 0.8. Errors are smaller than the symbols. The continues lines are drawn for 
guiding the eyes. 



and 



R f 



/• 



dA4A sech' 



7i / ,o a v 
— 2A 2 + — - — 
av \ 2 8 , 



(24) 



Inserting the parameters of the initial soliton (v o = 0.8, ao = 0.5) and then 
calculating the integral numerically (using Maple V) it gives the simple equa- 
tions R N ~ 1.31VQ 2 and R E — O.TIV^ 2 . This corresponds to a transmission 
coefficient for the mass 



Ti 



N 



1 - R 



N 



1.3T^ 2 



(25) 



and for the energy 



T; 



E 



1 - R K ~ 1 - 0.71K 2 



(26) 



The quanitative behavior of the transmission coefficients measured from our 
simulations is shown in Fig. 1. The height of the simulated potential step 
ranges from almost to about three times the kinetic energy of the soliton. 
The transmission coefficients fill monotonously the range from 1 to 0. We want 
to point out that the value of T E is always larger than T N . This is due to the 
fact that the kinetic energy is larger than the internal binding energy of the 
soliton (3i> > ao)- 

In Fig. 2 the result of perturbation theory (25) and (26) is compared with 
our numerically measured transmission coefficient T N for small parameters 
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Fig. 2. Comparison between the measured transmission coefficients T\r and Tg 
and the perturbation theory result Tjv = 1 — 1.31V^ 2 and = 1 — 0.71Vq 2 , which 
is expected to be valid for Vq/Kq <C 1/2. 
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Fig. 3. The transmission coefficients Te(B), and Tp(*) vs. Vo as in Fig. 1 

in a logarithmic scale. An exponential function is fitted to each set of the data 
(points) to indicate its behavior for large values of Vq/Kq. 



Vq. Due to (17) we expect that the perturbation theory is valid in a range 
Vq/Kq <C 1/2 in the present case. In fact, one sees a rather good agreement 
between the data and the first order perturbation theory for low step heights. 
The discrepancy for larger perturbations (Vq/Kq > 0.25) indicates that higher 
order perturbation theory would be necessary to predict the transmission of 
the soliton in this region more precisely. One has to bear in mind that the 
perturbation analysis assumes that the initial soliton does not break up and 
no soliton is reflected. 

In the opposite case, where the soliton is reflected at very high potential 
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steps (Vo Ko), we find an exponential decay of the transmission coefficients 
(see Fig. 3). All three transmission coefficients (Tjv, Tg, and Tp) show this 
exponential decay for V > K . The exponential functions plotted in Fig. 3 
are fitted by hand. This simple behavior strongly suggests that one should be 
able to derive it analytically. Up to now we did not yet succeed, but we are still 
working on this problem. Intuitively it is clear that the higher the potential is 
the less quasi-particles penetrate into region il. This amount of quasi-particles 
decreases exponentially with the penetrating depth of the soliton, because the 
tail of the soliton decays exponentially. Therefore the amount of quasi-particles 
"pushed" into and staying in region II decreases exponentially with increasing 
V . 

3.1.2 Potential ramp 

Let us now study potential ramps for positive values of Xq, i.e. potential ramps 
with finite slopes. Our data show unambiguously that the slope has a strong 
influence on the scattering process. In the equivalent particle method [36], 
where the soliton is treated like a classical Newtonian particle with mass iV 
and kinetic energy K , one would expect that the transmission coefficients 
remains constant when changing the length of the ramp. This would be due 
to the constant difference V — K . But it turns out that this simple picture 
is not true in our case. Bearing in mind that the transmission coefficients for 
a Gaussian wave packet of the linear Schrodinger equation being scattered 
at a potential ramp shows a dependence of the slope, too, this is not that 
remarkable. But the windows of reflection are. 

In figure 4 the transmission coefficient T N are plotted for different slopes and 
heights of the ramp. A similar figure with less complete data was already 
shown in Ref. [38]. In the following we will show mainly data for a initial 
velocity vq = 0.8 without loss of generality. We also have simulated different 
initial velocities vo and have found very similar results. 

Let us first shortly recapitulate the observations found in Ref. [38] concerning 
the global behavior of the transmission coefficient. After that we want to 
discuss the windows of reflection in a more detailed manner. 

The obvious stucture one observes is that the curves for Vq > K — 0.32 tend 
to if xq — > oo, while those for Vq < Kq tend to 1 in this limit. This indicates 
that the equivalent particle method is valid in the limit of very flat ramps. If 
the kinetic energy of the soliton is smaller than the potential height, it will 
be totally reflected; otherwise it will be totally transmitted. Simulations with 
still flatter ramps, e.g. xq = 100 (not shown here), show that the transmission 
coefficients are equal to 1 resp. within the error bars. 

Another property is the nonmonotonic behavior of the transmission coefficients 
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at steep ramps. All curves Kq < Vq show a global maximum at Xq ~ 3. 
If Vq is decreased below K , the absolute maximum is reached in the limit 
x — > oo (that is 1), but there is still a local maximum as long as V ^ 0.31. 
At some threshold this local maximum transforms into a shoulder. In Fig. 4 
this threshold is at V 0.31 for T N , while it happens at V ~ 0.3 for Tg 
[38]. For all values of Vq < 0.31 the position of this bump nearly stays at the 
same place. Simulations with different initial velocities (vq = 0.2, 0.6) likewise 
show such a maximum resp. shoulder at x ps 3. The fact that the maximum 
resp. the shoulder is located at a length of the ramp which is roughly equal to 
the width of the soliton, strongly suggests that this is a nonlinear resonance 
effect. A very similar behavior of the reflection coefficient was found by Yu.S. 
Kivshar et al. [37] for solitons scattered by a single impurity. 

Besides the global maximum, the curves for high potentials (Vo > 0.36) show a 
slight shoulder at 5 < x < 10. This shoulder gets more and more pronounced 
if Vo is decreased towards the magnitude of the kinetic energy (Vo = 0.33). At 
the same time its position is shifted towards flatter slopes (xq ~ 7 — > 12) and 
the shoulder becomes sharper. In case the height reaches K (V = 0.32) this 
shoulder is a dominating feature. 

Decreasing the potential further (below K , but still above some threshold 
Vo )C ), this shoulder continues to become sharper and more pronounced. In- 
deed, it develops into a peak and the transmission coefficients begin to grow 




Fig. 4. The transmission coefficient Tjv vs. xo for different potential heights and 

V = 0.8. 
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Fig. 5. The transmission coefficient Tn vs. xo as in Fig. 4, but displaying ad- 
ditional curves revealing windows of reflection (0.315 < Vq < 0.32) with higher 
resolution. 



with xq before they rapidly drop to join a common curve. This curve seem to 
be independent of Vq within the interval Vo )C < Vq < Kq. But the transmis- 
sion coefficients do not stay on this common curve as we increase xo further. 
We observe that the limit value for very flat slopes is not reached smoothly. 
Instead, the transmission coefficients jump back — instantaneously within our 
resolution — to values ~ 1. The values of x , where these jumps happen, seem 
to diverge when V — > K from below. 

As V is decreased further, we reach a value V 0jC = 0.3155 ± 0.0005 where the 
window of small transmission coefficients disappears. As Vo jC is reached from 
above, the locations of the downward jump and of the upward jump both 
tend towards xo, c ~ 14.7. At this point, both jumps seem to be instrumentally 
sharp. For V > Vo jC only the second jump retains this sharpness. For V < V 0)C , 
no trace of this singular behavior is left, and the transmission coefficients 
rise monotonically towards their asymptotic value 1 as xo is increased. An 
enlargement of the relevant part of Fig. 4, together with some more data, is 
shown in Fig. 5. 

These unexpected and surprising windows of small transmission coefficients 
were first seen in [38]. This effect is not expected within the picture of the 
equivalent particle method or within the linear Schrodinger equation. Super- 
ficially similar effects have been observed within the SG Eq. [39] and the 
4 -model [40] . But there are some fundamental differences to the actual work. 
These authors have studied the effect of a local and attractive impurity which 
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Fig. 6. Time evolution of the local maxima of |^| 2 for ten different slopes of the 
potential ranging from xq = 12.2 to xq = 16.0. The initial velocity is vq = 0.8 
and the height of the potential is Vq = 0.3180 = 0.99375i^o- Only the maxima 
|^| 2 > 1/3000 are plotted, and only the most important maxima are labeled. The 
calculation was done on a lattice with 8192 sites, discretization width Ax « 0.2 and 
integration step At = 0.005. 



is such that the soliton can be captured, and the energy is stored in an internal 
mode. This leads to sharp resonance frequencies. Another difference is that the 
soliton is unbreakable in [39,40]. Thus, in contrast to our system, the soliton 
is either reflected or transmitted or captured (in a state oscillating around the 
impurity) and does not emit radiation. The authors of [39,40] observe that the 
soliton may be reflected by the attractive impurity when its velocity lies in 
one of several windows. Within our resolution we find only one window, where 
the soliton is reflected instead of being transmitted. There are no indications 
for further windows of reflection. 

To understand this phenomenon in more detail, we first try to get some more 
insight by looking at the evolution of the wave function. In Fig. 6 we show the 
local maxima of |\l/| 2 for Vq = 0.3180 = 0.99375i^o- I n this diagram one can 
directly compare the effect of different slopes on the movement of the soliton. 
Let us begin with the steep ramp x — 12.2 (label 1 in the figure). Here 
the soliton lowers its speed while it climbs the potential ramp. At the same 
time radiation is emitted in the forward and backward direction. The biggest 
amount of radiation is emitted at the lower part of the ramp and immediately 
forms a new soliton escaping in the backward direction. The slowing down 
process continues until t ~ 350 and is accompanied with loss of radiation in 
the backward direction. After that the soliton moves with constant velocity 
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in the forward direction. But its mass has shrunk to roughly 60% of its initial 
mass. 

Flattening the slope of the ramp, the velocity of the transmitted soliton re- 
duces continuously. Whereas its mass decreases only slightly. If the slope is 
Xo = 12.8 the soliton interacts with the ramp until t ~ 1000. This corresponds 
to a distance x ~ 18 between the maximum of the soliton and the upper edge 
of the ramp. The same holds for the steeper slopes described above. 

Flattening the slope a tiny bit more (xq = 12.9) the soliton nearly stops after 
penetrating in to region II. It takes very long (t ~ 1400) until the soliton 
decides in what direction to move on. One consequence is that more radiation 
is emitted into region / while the soliton interacts with the inhomogeneity. 
The mass shrinks to roughly 40%, and it takes until t « 2700 that the soliton 
accelerates to full speed (not shown in Fig. 6). 

This behavior changes if the slope is again only slightly flatter. For x = 13.0 
the soliton moves with constant positive velocity already at t ~ 1100. But its 
loss of mass due to radiation is much higher. The transmitted soliton has only 
roughly 33% of its initial mass. This loss creates more than one new soliton 
and nonsolitary waves escaping from the ramp in the backward direction. It 
is hard to distinguish between solitons and nonsolitary waves here, because 
the amplitudes are small and furthermore the amplitude oscillates with a 
very low frequency, which makes it hard to identify them as solitonsQ One 
should notice, that nonsolitary waves are created near the lower edge of the 
ramp at least until t ~ 900. The other waves are not plotted, because their 
amplitude is less than 1/3000. As one can see from Fig. 7 the transmitted 
soliton continuously loses mass while it interacts with the potential ramp. This 
happens until it leaves the upper edge of the ramp and moves with constant 
speed again. 

If one looks closer at the traces of the transmitted soliton in Fig. 6 one observes 
small bends at times t w 200, 350, 500 and 650. It looks as if the soliton's 
maxima are oscillating around a true trace. This frequency corresponds exactly 
to the oscillation period of the amplitude of those solitons. More precisely this 
means that the apparent bend in the trace occurs when the transmitted soliton 
is maximally sharp. This can be explained in the following way: the soliton 
is stretched when it climbs up the ramp, which causes the amplitude to be 
reduced. Thereby it becomes asymmetric as one can see in Fig. 2 of Ref. [38] 
and its tail stays on the ramp for a long time. During that time the pulsating 
soliton still interacts with the ramp, which causes the oscillations seen in Fig. 6. 

While decreasing the slope of the ramp more and more the mass of the trans- 

1 The reader is refered to section 4.1 where the amplified solitons are discussed in 
more detail. 
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Fig. 7. Time evolution of the transmission coefficient T/v of four representa- 
tive values of xo (12.2, 12.8, 13.0 and 14.0). The height of the potential is 
V = 0.3180 = 0.99375Ko as in Fig. 6 



mitted soliton resp. solitary wave drops rapidly but continuously to ~ 10% 
(see Fig. 7). This is the reason why the speed of the soliton increases at the 
same time. For slopes 13.5 < x < 15.0 the movement of the penetrating soli- 
ton is not monotonical. After it has come to a stop at t ~ 200 it moves back 
before it stops again and escapes in the forward direction. E.g. let us look at 
the situation with x = 14.0: at t ps 200 the soliton becomes wider, which can 
be concluded by comparing the oscillations of the other traces shown in Fig. 6. 
Due to this broadening its maximum is shifted towards region I. During this 
process some of its mass is pushed down the ramp into region /. That can be 
concluded from Fig. 7 where the curve shows a drop at T ps 200, after it had 
seemed to converge to a value larger than that for x = 12.8. After some part 
of the soliton has separated from the bulk it accelarates away from the ramp 
while it contracts again. 

For xo = 16.0 the soliton penetrates into region 77, but is reflected at the top 
edge of the ramp. This is consistent with the collective-coordinate description 
(see Eq.(ll)). The reflected soliton consists of roughly 96% of the initial mass. 
This scenario holds for flatter slopes, too, until Xo = 19.1. There the soliton 
is transmitted with roughly 95% of its initial mass and only few radiation is 
emitted. 

For heights V close to V 0jC the first jump becomes sharper and finally as sharp 
as the second one. 

As already pointed out in our previous work [38] the second upward jump can 
be explained easily: If there were no radiation (as is the case for x — > oo) the 
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energy of the soliton would be conserved, and the threshold for transmission 
would be strictly at Vo = K . But for finite Xq, the radiation implies that the 
soliton has lost energy when it reaches the upper end of the ramp, and it needs 
K > V to be transmitted. Thus there is a range K — S(x ) <V <K where 
the soliton is reflected. Since the soliton moves very slowly near the upper egde 
of the ramp for Vo ~ K$ — 5(xq), the dependence of the transmission coefficients 
on the amount of radiation (and thus also on xq) is extremely sharp. To be 
more quantitative, let us consider V = 0.3190. Here, the second jump (from 
nearly total reflection to nearly total transmission) occurs at x = 22.1 = x c . 
At this slope, T E jumps from 0.0093 to 0.9837, while T N jumps from 0.0074 
to 0.9797. The fact that the soliton is coming to a near stop for xq = x c + e 
means that all its energy is potential energy, E' = N'(Vq — a' 2 /6), and the 
energy difference has gone into radiation. Primes denote here observables at 
the moment when the soliton stops. If we neglect the forward radiation (which 
seems to be indeed small according to Figs. 6), we have N' = T N , E' = T E E 
and a' = (T N /2) 2 , and thus T E = T N (V - T^/2A)/E = 0.9821 for x = x c . 
Comparing this with the measured value 0.9837 we see good agreement. 

Within the window of reflection, the time it takes until the point of return 
is reached is not the same for different slopes. This time becomes larger and 
larger the nearer the slope is to the border of the window. Indeed, that is 
true for both borders. The reason for this is that the radiative loss of energy 
is equal to the difference K — Vo. Therefore the velocity of the temporarily 
transmitted soliton is very close to 0. Thus, the soliton stays during a long 
time in the vicinity of the egde of the ramp and can interact with it. This 
means that already tiny amounts of radiation can cause the soliton to turn 
back into region /. The same argument holds when decreasing the difference 
Kq — Vo- The smaller it is, the longer it takes until the soliton has decided 
whether to pass the potential ramp, or to fall back into region /. 

In order to exclude that the observed feature is due to the non-differentiability 
of the chosen potential, we have simulated smoother potentials, too. We found 
that the values of the transmission coefficients change only slightly if one uses 
instead of potential (2) a more smooth sigmoid. Even the shape of the windows 
with small transmission are still seen and they stay at the same positions. 
However, if the potential becomes too smooth (e.g. such as functions like 
tanh(srr) with s < 1), the soliton behaves as one expects for the equivalent 
particle method and only very little radiation is emitted. 



3.2 Transmitted momentum 

Another observable we can look at are the momenta of the reflected and trans- 
mitted wave. Although momentum is not conserved for nonconstant potentials, 
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it is worth to have a closer look at it. If the wave function is only composed of 
a single soliton, its momentum is directly related to its velocity via Eq.(7b). 
In order to measure the momentum of one soliton within a system of several 
nonoverlapping solitons, it is necessary to restict the integration interval (in 
Eq.(6)) around its position. If we neglect the momentum of the radiation, and 
if at most a single soliton is emitted in the forward and backward direction, 
it is sufficient to partition the a;-axis into regions I and 77. That is what we 
have done. 

As mentioned before, according to the equivalent particle method the trans- 
mitted momentum should not depend on the slope of the ramp. But what we 
observe is different, as already pointed out in the previous section. The global 
behavior of the transmission coefficient is the same as it is for the transmit- 
ted mass and energy. In Fig. 8 we see a (local) maximum at steeper ramps 
(xo ~ 2 — 3), which confirms the claim that this is a nonlinear resonance effect. 
In the limit Xo — > oo the momentum tends to the value one expects for the 
equivalent particle method. There one assumes that if the soliton is transmit- 
ted, it is so without loss of radiation. Due to the conservation of energy the 
reduced velocity is given by vn jOUt in the cases where V < K Q . Otherwise the 
soliton is reflected and the momentum of the transmitted solitary wave tends 
to zero. Our results agree with this expected behavior for very flat ramps 
(x > 25). 
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Fig. 8. The transmission coefficient Tp vs. xq for different potential heights as in 
Fig. 4. 
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Fig. 9. Enlargement of Fig. 8 on a logarithmic scale. Only transmission coefficients 
for potential heights 0.3150 < Vq < 0.33 are plotted. 



Between the (local) maximum and the limit region (x — > oo) there exists a 
minimum whose position depends on the height of the potential. For lower 
heights it is located at steeper ramps (xq ~ 8). While raising the potential, 
the minimum gets more pronounced and is shifted to flatter ramps (x ps 14). 
When the ramp is nearly as high as the kinetic energy of the soliton (Vo < 
0.3150), the values of the transmitted momentum at the minimum is very 
small and the minimum becomes a dominating feature. 

Raising the potential above the threshold Vq jC (~ 0.3155; see subsection 3.1.2) 
the situation is different. For V G [Vo, c , Ko] the momentum of the transmitted 
soliton has very sharp minima at both ends of the window of reflection (see 
Fig. 9). In the interior of the window it is non-zero and decreases exponentially 
with Xq. The minimum at the upper end of the window tends to zero, and the 
same holds for the lower end if Vo is not too close to K$. But for Vo ~ Kq 
(Vo > 0.3185) when also the lower end of the window of reflection in Fig. 5 is 
slightly smoothened, T P has a smooth dip of finite depth. The larger V is, the 
smoother the behavior gets. Finally, it vanishes at a height V > 0.325. Even 
the curve for V = K = 0.32 shows a smooth dip at x 12.7. The cusp at 
the second edge remains sharp for values 0.3150 < Vo < 0.3200. Whereas the 
first cusp shifts only little from xq = 14.1 to xq = 12.7 as Vo — > Kq, the second 
cusp shifts rapidly from x = 14.9 to x > 22. For potentials higher than the 
kinetic energy we do not find any jump from T P rs to T P « 1 at all. 
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Outside the window but still in its vicinity, the mass of the transmitted soliton 
stays nearly constant. This leads to the conclusion that the velocity of the 
transmitted soliton drops to almost zero at the edge of the window. This 
agrees with the observation made in Fig. 6. The dependence of T P inside the 
window is dominated by the mass of the transmitted solitary wave, because 
its velocity changes only slightly. It is conspicuous that all curves within the 
window match the exponentially decaying curve for Vq = 0.33. This confirms 
the conclusion that those transmitted solitary waves are universal in the sense 
that its mass, energy and momentum are independent of the height of the 
potential. 



4 Shape of the reflected and transmitted soliton 

4-1 Description by an amplified soliton ansatz 

After the soliton hits the potential ramp, we observed in general that it splits 
up and creates a transmitted and reflected soliton, which are both accompa- 
nied by some radiation (see Fig. 10). This is why we find typically more than 
a single maximum of |\l/(a;)|. Moreover, the heights of these maxima in general 
are not constant in time, as one would expect for a pure soliton state. Instead, 
the amplitude of the soliton shows in general marked oscillations (e.g. Fig. 11) 
which are damped in all cases. In spite of this, the solitons moves with practi- 
cally constant velocity (see Fig. 10; the small bent lines belong to the maxima 
of nonsolitary waves). 

Damped amplitude oscillations have be observed in connection with 'amplified' 
solitons. These have been studied analytically and numerically in [41,42,34]. 
These authors have shown that a soliton whose amplitude is multiplied by 
some factor relaxes by emitting radiation and by showing damped oscillations 
of its amplitude. This can be understood by interpreting the initial state as a 
superposition of a soliton with radiation. 

As an example we just show one typical scenario to illustrate the movement 
and the amplitude behavior of the solitons during the scattering process. Other 
scenarios are shown in our previous work [32]. Figures 10 and 11 depict the 
case where the potential is a step function (x = 0) and the kinetic energy 
(K = 0.32) is larger than its height V = 0.3. In the equivalent particle 
method one would expect the soliton to move into region II and to propagate 
there with a reduced speed vn,out- But our simulation shows that it breaks up 
into two solitons with roughly equal heights and with velocities v = —0.588 
and w = 0.395. About half of the mass and three quarters of the energy are 
transmitted (T N = 0.527, T E = 0.712). 
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The loss due to direct radiation is very small. If one inserts the above numbers 
into Eq.(7c) one finds perfect agreement (discrepancies are < 1%) [32]. This 
indicates that radiation in form of nonsolitary waves is small in spite of the 
wiggles seen in Fig. 11. 

Assume now that the wave functions after the interaction can be described in 
both regions separately by amplified solitons 

*(x',t> = 0)= a + 6 ,^ x \ (27) 
cosh(arr) 
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Fig. 10. Time evolution of local maxima of 1^1 2 for a soliton with incident ve- 
locity vq = 0.8 which is scattered at a potential step with xq = and height 
Vq = 0.3 = 0.937^0- The calculation was done on a lattice with 4096 sites, dis- 
cretization width Ax = 0.2 and integration step At = 0.005. The latter parameters 
are the same for the next figures. 
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Fig. 11. Time evolution of the height of the maxima shown in Fig. 10. The highest 
curve belongs to the transmitted soliton, and the second highest to the reflected one. 
The other maxima presumably are due to the superposition of nonsolitary waves. 
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where x' = x — x and t' — t — t. The parameters x, t specify the location and 
the time of the interaction and v the new velocity. This function is a soliton 
whose amplitude is simply multiplied by a factor of 7 = 1 + f ■ The long- 
time asymptotic solution consists just of a soliton. For finite times it shows 
oscillations [41,42]. They are damped, and the soliton's amplitude tends to 

k = a + 25 . (28) 



The mass of the amplified soliton (27) at time t' — is given by 

(a + 5f 



:= [dx'\^(x',t > = 0)\ 2 = 2^-^- (29) 

J CI 



— oc 



and its maximum is |^(0, 0)| = a + 5. Satsuma and Yajima [41] proved that 
the nonsoliton part of Eq. (27) behaves like an ordinary wave packet and thus 
decays as t~ l l 2 for t — > 00. The mass of the asymptotic soliton is given by 

00 

:= Jdx' \V(x',t' = oo)| 2 = 2(a + 2<5) , (30) 

—00 



the difference with Eq. (29) being due to the unobservable radiation. 

Now we will use Eqs. (29) and (30) to evaluate the parameters a and S for the 
scattered wave functions. Then we will compare the scattering data with a 
simulation using as initial condition a wave function (27) with the calculated 
parameters. 

The parameters x, t, v have to be chosen in such way that the amplified soliton 
fits to the scattering data. The velocity v of the scattered soliton can easily 
be extracted from a plot that shows the time evolution of its maxima (see e.g. 
Fig. 10). However, the values of x and t are not fixed a priori, because they 
depend on the length of the ramp. 

In the case of the potential step (xq = 0) it is natural to place the interaction 
point simply at x — 0, t — —x s/vq, where x g < is the position of the 
soliton at time t — (see Eq. (16)). But if xq > the coordinates of the 
interaction point are different from those. 

To obtain the parameters a and 5 of the amplified soliton, the equations 

iV = 2(a + 25 + 5 2 /a) (31a) 
iV^ = 2(a + 25) (31b) 



23 



have to be solved. Here, iVo is simply either the transmission or the reflection 
coefficient, T/v or R^. The value of is taken from plots which show the 
time evolution of the amplitude. From them it is easy to estimate k = a + 25, 
but the determination of a and 8 themselves is less straightforward. 

Solving Eqs. (31) we get two sets of solutions: 

ai,2 = N -^± Jn (N - Nn) , (32a) 



Noo - N =f JN (N - iV^) 
81,2 = — V 2 1 , (32b) 

corresponding to initial amplitudes which are larger (Si < 0) resp. smaller 
(8 2 > 0) than that of the asymptotic soliton. Also, the period of their oscilla- 
tions are slightly different. 

Unfortunately, those relations do not help to decide which solution, (a l5 Si) or 
(a 2 , S 2 ), describes the scattered soliton. This can be found out by direct com- 
parison, i.e. numerically, because there does not exist an analytical solution 
of the NLSE with initial condition (27). As an example, we plotted the two 
solutions together with the simulations in Figs. 12 and 13, for v = 0.8 and 
V = 0.3 

As one can clearly see, only the amplified soliton with ai = 0.4894, Si = 
—0.1304 is in good agreement with the data. Both amplified solitons show 
oscillations with almost the same period as the transmitted one. But the am- 
plitude fits much better for the first solution, and the behavior around the in- 
teraction point is much better represented. Also, the side maxima are not that 
close to the main maximum as it is for the curve with a 2 = 0.1067, 82 = 0.0609. 
All this is an stong indication that the transmitted soliton can be described 
by an amplified soliton ansatz with solution ai,Si. 

Even better fits are obtained in Figs. 12 and 13 if the parameter a is not 
calculated from Eqs. (32), but is simply chosen as the value a = = 0.5 of 
the unperturbed soliton. Within errorbars, this is consistent with the above 
value. The corresponding 8 is —0.13575. 

For this case (t>o = 0.8, Vq = 0.3, xq = 0) the interaction of the incoming soliton 
with the potential can be interpreted in the following manner: the soliton with 
amplitude 0.5 interacts with the potential inhomogeneity (xq = 0). The main 
part (T E Rs 0.7) is transmitted, and the remaining part is reflected. Most 
of the radiation is emitted in the backward direction. Due to the birth of a 
second soliton the amplitude of the initial soliton is reduced by 5. This new 
state of the wave function is not a pure soliton any more. Thus, the amplitude 
decreases further and begins to oscillate in the manner seen in the figures. The 
reflected wave also forms such an amplified soliton. In addition, radiation is 
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Fig. 12. Comparison between the amplitudes of the two simulated solu- 
tions ai = 0.4894 ± 0.0400, ft = -0.1304 ± 0.0200, a 2 = 0.1067 ± 0.0150, 
82 = 0.0609 ± 0.0070, and the transmitted soliton, Fig. 11) with the parame- 
ters v = 0.8, x = 0, V = 0.3. 
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Fig. 13. Comparison between the position of the local maxima of the two simu- 
lated solutions and the transmitted soliton as in Fig. 12. 

created as seen from the various decreasing and non-oscillating side maxima. 

For natter slopes of the ramp, amplified solitons fit well for times t' larger than 
one oscillation period of the amplitude. The early stages of the interaction can 
not be reproduced with this ansatz. In our simulations we observed that the 
amplitude oscillation is slightly compressed horizontally while the soliton is 
located in the rising part of the potential. This might be because of the transfer 
of energy into the radiation. 

The amplitude oscillations for the reflected soliton are less easy to compare. 
This comparison is possible for very steep ramps and rather high potentials. In 
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Fig. 14. |^| 2 of the simulated amplified soliton with m = 0.6461 ± 0.0400, 
5\ = —0.1006 ± 0.0100 in comparison with the reflected soliton with the param- 
eters Do = 0.8, xo = 0, Vq = 0.5. 
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Fig. 15. Comparison between the position of the local maxima of the simulated 
solution and the reflected soliton as in Fig. 14. 

those cases the reflected soliton has a larger amplitude than the nonsolitary 
wave emitted during the interaction. Thus, effects due to interferences are 
rather small. If their amplitudes are of the same order, they strongly interfere 
and deform the shape of the amplified soliton. A way to avoid that problem is 
to wait very long until the radiation either has passed through the soliton or 
is dispersed. But this would require very long simulations. Another problem 
is that in general the oscillation periods of the reflected soliton get very large 
which requires a very long simulation time, too. 

Another comparison between data and theory is shown in Figs. 14 and 15 
for the situation v = 0.8, V = 0.5, and x = 0. The simulation and the 
amplified soliton show only small discrepancies, mainly during the first three 
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oscillation periods. The interference is also the reason why the side-maxima 
of the radiation emitted by the amplified soliton are not reproduced correctly. 
The rest is quantitatively very well reproduced. 

In summary, each of the solitary waves after the scattering can be described as 
an amplified soliton, except for a very short time during the collision process. 
It is described by means of only two parameters: the reflected resp. transmitted 
mass and the height of the soliton in the limit t — > oo. Since these states are 
not solitons at finite times, they emit radiation. This happens continuously 
and the amount radiated until a time t is proportional to Vt -1 [41]. Indeed, 
we see that the amplitude of the linear waves does not decay exponentially 
[42]. In the limit t — > oo, the solitary waves converge to solitons with defined 
shapes. The total amount of radiation lost is simply 

N -N oo = 2-. (33) 
a 

But this means that the overall amount of radiation produced during the whole 
scattering process is in general much bigger than one might have supposed at 
first sight. For steeper ramps more than 10% of the mass and even slightly more 
of the energy is eventually radiated by the two amplified solitons leaving the 
ramp into two directions. Summing up all effects we find that in cases where 
the potential is of the same order as the soliton's kinetic energy and the slope 
is steep, about one fourth of the initial mass and energy is radiated. Otherwise 
the loss is smaller. The percentage of radiation decreases systematically with 
decreasing slope. In the limit of very flat ramps (x — > oo) it tends to 0. 

4-2 Oscillation period of the solitons amplitude 

Kath and Smyth [34] derived two formulas which estimates the period of 
the amplitude oscillations for an amplified soliton with shape (27). For this, 
they used an alternative method to describe the time evolution of the soliton, 
namely a Lagrangian approach [43] . This leads to four approximate differential 
equations. These coupled ODE's allow to predict the time evolution of the 
soliton's parameters including the "mass loss" due to dispersive radiation. 
Their results are in agreement with the full numerical solution of the NLSE. 

As we have shown above, the assumption that the solitary wave after the 
scattering is described by an amplified soliton with suitable values for a and 
5, seems to be valid in a wide range of external perturbations. Thus we can 
apply their formulas to our data to see how far they agree. 

For small perturbations, it was shown in [34] that the period of the amplitude 
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Fig. 16. Oscillation periods of the transmitted soliton (data), compared with the 
ones calculated using Eq.(34). Small values of xo correspond to large \S\ and large 
xo to small \5\. The initial parameters of the system are Vq = 0.30, vq = 0.8 and 
ao = 0.5. The errors on the directly measured data are not shown, but are roughly 
half as big as those plotted for Eq.(34). 

oscillations is the same as that of the phase oscillations of the pure soliton at 
time t — oo: 



For finite perturbations (\S\ > 0), the same authors found (using our notation) 



Notice that this does not agree with Eq.(34) for \S\ — > 0, as the authors of [34] 
already pointed out. 

Measured oscillation periods are plotted in Figs. 16 and 17 against x for 
two different velocities. These are comparisons with the predicted values of 
Eq.(34), where we have substituted k = a + 25. The values predicted by 
Eq.(35) do not fit the data in a sufficient manner. It predicts smaller values 
for x > 16 (v = 0.8) resp. x > 10 (v = 0.6) and much larger values for 
x — > 0. But the simple Eq.(34) shows a surprisingly good agreement for all 
data points in both cases. 

We could conclude from this that the disturbance of the soliton's amplitude 
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Fig. 17. The same as in Fig. 16 but for different initial parameters (Vq = 0.17, 
vq = 0.6, ao = 0.5). 

we are dealing with seems to be not as large as expected. Namely, in the sense 
that Eq.(34) — meant to be valid only for very small disturbances — predicts 
the correct time periods. This is astonishing, because we found an reduction 
of the soliton's amplitude up to 25%! This would mean that the valitiy range 
of Eq.(34) is much larger than assumed in [34]. 



5 Summary and conclusions 

In this paper we found that solitons break up in general when hitting a poten- 
tial threshold. If the slope is very steep and the potential height is either very 
low or very high compared to the soliton's kinetic energy (Kq ^> Vq or Kq <C 
Vo), the soliton acts like a Newtonian particle. The same behavior is seen if 
the slope is rather flat (Vq/xq <C 1) for most values of Vq. This includes the 
fact that in those cases only very little energy and mass is transfered into 
radiation. In such cases only one soliton with nearly constant height exists 
during the whole simulation time. The amplitude of this solitons shows small 
oscillations with a period, that is the period of the phase oscillation of the 
undisturbed soliton. 

The complexity of the outgoing state depends on the parameters of the po- 
tential and of the soliton, but most frequently the soliton breaks into two 
solitons of different mass. We have shown that these solitons can be described 
by amplified solitons whose parameters a and 5 can be evaluated from the 
data. Additionally, the period of the amplitude oscillation can easily be calcu- 
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lated. The scattered solitons are accompanied by radiation, which is emitted 
continuously and whose amplitude decays proportional to 1/ \/t in time for 
large times. The fact that the radiation decays only with a power law raises 
some problems with the transmission of information via optical fibers [42]. 

In some region of the (xq, Vo)-plane the soliton is almost fully reflected even 
if its kinetic energy is slightly larger than the potential. This is due to the 
emission of radiation. The energy of the soliton is reduced thereby, whence 
it cannot overcome the ramp and is reflected. The associated jumps in the 
transmission coefficients are rather sharp and the slope where this happens 
depends on the hight of the potential. The window in slopes where the soliton is 
reflected shrinks if the height of the potential is lowered. At a critical value Vo )C 
this feature disappears and the transmission coefficients rise monotonically. As 
the potential height tends to the value of the soliton's kinetic energy, the right 
border of the window tends to +00. 

In this paper we have mostly studied the NLSE with an external potential 
that has the form of a linear ramp with variable hight and slope. We have ap- 
plied an optimized fourth-order symplectic integrator in our simulations. The 
integrator is optimized in the sense that it takes into account the symplectic 
structure of the NLSE and the fact that the kinetic energy is bilinear in ty x . 
It has a very good long-time stability. 

It is very easy to apply this algorithm to NLSEs with other external poten- 
tials, such as one or more impurities, a series of step functions or any smooth 
function. In the same way other Hamiltonian disturbances can be treated, 
e.g. terms like |\l/| 2p \l/. In order to study the collision of two solitons (a fast 
one and a slow one) or of three solitons (a fast one crossing two interacting 
slow ones) [44,45], an integrator with good long-time stability is needed. We 
believe that the integrator used in the present work can be very useful in 
further simulations of partial differential equations of such a type. This nu- 
merical method can also be easily and staightforward applied to the NLSE 
with a varying nonlinear term or to higher dimensional and coupled NLSEs. 
Thus the vector NLSE, which describes e.g. the propagation of light in a bire- 
fringent optical fiber, can be simulated with such an algorithm. Also soliton 
switching and splitting can be simulated by using this algorithm. 
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Appendix 

A Symplectic integration 

Here we briefly describe the integration method which we have used to simulate 
the NLSE. 

As we have already mentioned, the NLSE is a Hamiltonian system. Thus, 
it is natural to apply to it integration routines which were developed dur- 
ing the recent years and whose main characteristic is that they preserve the 
Hamiltonian structure [46-48]. The latter is not true e.g. for standard meth- 
ods as e.g. Runge-Kutta or predictor-corrector. Such 'symplectic' integrators 
(the simplest of which is the well known Verlet or 'leap frog' algorithm) have 
been applied already to the linear [49-52] and nonlinear [32,53-56] Schrodinger 
equations. 

The most popular algorithms of this type are split-operator methods. They 
depend on the Hamiltonian being a sum of two terms A and B, each of which 
can be integrated explicitly. Then one uses the Baker-Campbell-Hausdorff 
theorem to approximate e l ( A + B ) t by a product of factors e lolkAt and e^ 5 *, 
where and /3k are real numbers satisfying among others J2k a k = J2k flk = 1- 
The error is given by higher order commutators of A and B. In particular, we 
apply a fourth-order method introduced by McLachlan and Atela [57] which 
is applicable if one of the third order commutators vanished identically. We 
found that this method should be applicable to our problem, and that it is 
indeed numerically very precise, indicating that the McLachlan- Atela method 
is the method of choice for a wide class of problems. 

The NLSE is a classical Hamiltonian system with Poisson bracket 



and Hamiltonian H = E. This implies in particular that it can be written as 



where the linear ('Liouville') operator Ti is defined as Ti • = { ■ ,H}. Split- 
operator methods can be applied by splitting Ti = T + V, where T and V are 



{V*(x),V(y)} = i6(x-y) 



(A.l) 



* = = m, 



(A.2) 



31 



the Liouvilleans corresponding to ~ / dx^ 1 !/! 2 and / dx(— ||\l/| 4 + V|\I/| 2 ), 

T^ = -<9 2 #, = i{\m\ 2 ^ - Vtf) . (A.3) 

In [57], a fourth-order algorithm was introduced which minimizes the neglected 
fifth order terms in the Baker- Campbell-Hausdorff formula for Hamiltonians 
for which 

[[[T,V],V],V] = 0. (A.4) 

This applies obviously to Hamiltonians with T = |(p, M~ l p), V = V(x), with 
M a constant mass matrix and {qi,Pk} = $ik-, since in the Baker-Campbell- 
Hausdorff formula each commutator with V acts as a derivative operator on 
any function of p. In [51] it was shown that this algorithm can also be applied 
to the linear Schrodinger equation where it gives better performance than 
the general fourth-order algorithm [46] which does not take into account this 
special structure. In [32] it was shown that eq. (A.4) holds also for T and V 
defined in eq. (A.3) and the McLachlan-Atela method can thus also be used 
for the NLSE. 

The coefficients and (3k for the McLachlan-Atela method are listed in 
[57,51]. Our implementation involves a spatial grid with Fourier transforma- 
tion after each half step [51]. 

Since T and V both conserve the mass exactly, N should be conserved up to 
round-off errors. This was checked numerically, relative errors typically were of 
order 10~ n . Energy is not conserved exactly, and its error was ~ 10~ 7 after an 
evolution time t = 500 with an integration step At = 0.0025. The precise value 
of the error depends of course on the parameters of the soliton and on xq. We 
checked carefully that our final results were independent of the time step and 
of the spatial discretization Ax. It was checked that the algorithm is indeed 
fourth-order, and is more precise than the general fourth-order symplectic [46] 
and the leap-frog (second-order symplectic) algorithms. 

The derivatives of \I> and \&*, which have to be calculated numerically, were 
computed in Fourier space, as this is much more precise than taking finite 
differences in x-space. Using the latter, the relative error of the energy con- 
servation would have been of the order 10 -4 instead of 10 -7 . 



2 As an alternative to the Fourier transform one can perform these steps also in 
x-space [58]. Instead of a FFT, one has to perform a convolution with the Greens- 
function of T in each time step. Practically, both methods are roughly equally fast. 
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